  	setwd("C:/Users/Jef/Desktop/Applied Statistics")
#  	setwd("D:/Mijn Documenten/applied statistics - project")
	library(qcc)

# voor standaard normale variabelen:
n		<-	100
d1N		<-	qnorm(1-(-0.7+0.5*log(n))/n)
d2N		<-	qnorm(1-(5/(n*n^(1/2))))
Norm <- 0
Param <- 0
Minimal <- 0

for(i in 1:10000){
  #trek een sample van 100 standaard normaal verdeelde rv
  NormSample <- rnorm(n)
  #orden deze op grootte : X_(1)<...<X_(n) 
  SortNormSample <- sort(NormSample)
  #bereken mean en S (standard deviatie van deze 100 stochasten)
  Mean <- mean(SortNormSample)
  StDev <- sd(SortNormSample)
  #als d1N<=(X_(100)-mean)/S<=d2N dan is deze normaal verdeeld
  if(d1N<=(SortNormSample[100]-Mean)/StDev && (SortNormSample[100]-Mean)/StDev<= d2N){
    Norm <- Norm + 1
  }
  else{
    Gamma <- 1.1218 * log((SortNormSample[96]-Mean)/(SortNormSample[76]-Mean))-1
    GammaFunctie <- pi^(1/4)*2^(-(1+Gamma)/2)*gamma((Gamma+3/2)^(-1/2))
    d1P <- GammaFunctie*(qnorm(1-(-0.2+0.5*log(n))/n))^(1+Gamma)
    d2P <- GammaFunctie*(qnorm(1-3/(n*n^(1/2))))^(1+Gamma)
    if(d1P<=(SortNormSample[100]-Mean)/StDev&& (SortNormSample[100]-Mean)/StDev <= d2P){
      Param <- Param + 1
    }
    else {
      Minimal <- Minimal + 1
    }
  }  
}
Norm
Param
Minimal

# voor exponentieel verdeelde variabelen mean 1:
n		<-	100
d1N		<-	qnorm(1-(-0.7+0.5*log(n))/n)
d2N		<-	qnorm(1-(5/(n*n^(1/2))))
Norm <- 0
Param <- 0
Minimal <- 0

for(i in 1:10000){
  #trek een sample van 100 exponentieel verdeelde rv's met mean 1
  NormSample <- rexp(n)
  #orden deze op grootte : X_(1)<...<X_(n) 
  SortNormSample <- sort(NormSample)
  #bereken mean en S (standard deviatie van deze 100 stochasten)
  Mean <- mean(SortNormSample)
  StDev <- sd(SortNormSample)
  #als d1N<=(X_(100)-mean)/S<=d2N dan is deze normaal verdeeld
  if(d1N<=(SortNormSample[100]-Mean)/StDev && (SortNormSample[100]-Mean)/StDev<= d2N){
    Norm <- Norm + 1
  }
  else{
    Gamma <- 1.1218 * log((SortNormSample[96]-Mean)/(SortNormSample[76]-Mean))-1
    GammaFunctie <- pi^(1/4)*2^(-(1+Gamma)/2)*gamma((Gamma+3/2)^(-1/2))
    d1P <- GammaFunctie*(qnorm(1-(-0.2+0.5*log(n))/n))^(1+Gamma)
    d2P <- GammaFunctie*qnorm(1-3/(n*n^(1/2)))^(1+Gamma)
    if(d1P<=(SortNormSample[100]-Mean)/StDev&&(SortNormSample[100]-Mean)/StDev <= d2P){
      Param <- Param + 1
    }
    else {
      Minimal <- Minimal + 1
    }
  }  
}
Norm
Param
Minimal




 




